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We identify qualitative trends in the stacking sequence dependence of carrier-carrier interaction 
phenomena in multilayer graphene. Our theory is based on a new approach which explicitly exhibits 
the important role in interaction phenomena of the momentum-direction dependent intersite phases 
determined by the stacking sequence. Using this method, we calculate and compare the self-energies, 
density-density response functions, collective modes, and ground-state energies of several different 
few layer graphene systems. The influence of electron-electron interactions on important electronic 
properties can be understood in terms of competition between intraband exchange, interband ex¬ 
change and correlation contributions that vary systematically with stacking arrangement. 


Introduction — Multilayer graphene has attracted con¬ 
siderable attention recently because of exotic chiral fea¬ 
tures in its electronic structure and the possibility of fu¬ 
ture electronic device applications [D-Q- The band struc¬ 
ture of a multilayer system is qualitatively dependent on 
its stacking sequence, opening up the possibility of en¬ 
gineering electronic properties by selecting a desired ar¬ 
rangement. In this paper we use an approach in which 
momentum-direction dependent intersite phases deter¬ 
mined by the stacking sequence are explicitly exhibited 
to show that this qualitative dependence is inherited by 
carrier-carrier interaction phenomena. 

Because the number of 7r-bands in a multilayer 
graphene system is proportional to the number of layers, 
and because 7r-band wavefunctions are not isotropic in 
momentum space, accurate evaluation of physical quan¬ 
tities which require integrations over momentum space, 
for example quasiparticle energy spectra and density- 
density correlation functions, rapidly becomes more diffi¬ 
cult as layer number increases. To mitigate this problem 
and to make the relationship between stacking arrange¬ 
ment and interactions more transparent, we introduce a 
momentum-direction dependent unitary transformation 
which makes the single-particle Hamiltonian isotropic. In 
addition to making accurate many-electron perturbation 
theory calculations practical for multilayer stacks, this 
approach facilitates understanding of some qualitative 
trends in the stacking arrangement dependence of quasi¬ 
particle energy spectra, plasmon dispersion and damping, 
and carrier thermodynamic properties. 

Rotational transformation of multilayer Dirac Hamil¬ 
tonian — Our calculation is based on the minimal con¬ 
tinuum model for multilayer graphene which retains only 
a Dirac model for hopping within each layer and only 
nearest-neighbor interlayer hopping. Different stack¬ 
ing sequences are specified by different interlayer near¬ 
neighbor arrangements. The Hamiltonians for these 
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minimal models can be made isotropic by multiplying 
wavefunction components by stacking and momentum- 
direction dependent phase factors. To illustrate how this 
transformation works, we consider first the example of 
Bernal stacked bilayer graphene, in which one sublattice 
in the first layer (say IB) is a near neighbor of the oppo¬ 
site sublattice in the second layer (say 2A). The Hamilto¬ 
nian at finite wavevector k is then expressed in the (1A, 
IB, 2A, 2B) basis as 

f 0 hvke~ l(t ’ k 0 0 \ 

hvke l ^ k 0 t± 0 

H[(pk) - 0 t_ L 0 hvke~^ k ’ 

y 0 0 hvke l<l>k 0 J 

(1) 

where k = y//c 2 + fc 2 , <j) k = arctan(fc y /A; x ), v is the bare 
Dirac velocity, which is related to the nearest-neighbor 
intralayer hopping amplitude by t = 2hv/\/3a ~ 3 eV 
(a = 0.246 nm is the lattice constant), and t± ~ O.li 
is the nearest-neighbor interlayer hopping parameter. It 
is easy to see that the eigenvalues of this Hamiltonian 
are independent of <t> k and that all eigenvalues satisfy 

'h(0fc) = (ciA,Ci B e^ fc ,C2Ae I ' ?ife ,C2Be 2 ^ fc ) t = U(c/>k)' k(0), 

where the {cj} depend on k only and can be obtained by 
diagonalizing H at (f>k = 0. The locking between inter¬ 
site phases and momentum direction in these spinors is 
reminiscent of the properties of spinors in chiral systems 
and will be referred to below as sublattice pseudospin 
chirality. The unitary operator U{(jf) is a diagonal ma¬ 
trix whose diagonal components (1, e*^, e 1 ^, e 21 ^) are de¬ 
termined by the bilayer stacking. The phase difference 
e 1 ^ between the 1A and IB components of the wavefunc¬ 
tion, and between the 2A and 2B components, comes 
from the monolayer-like intralayer coupling, whereas the 
zero phase difference between the IB and 2A components 
comes from the momentum-independent interlayer cou¬ 
pling. If we know the wavefunction at a specific angle, we 
can easily obtain the wavefunction at an arbitrary angle 
by attaching site-dependent phase factors determined by 
the stacking sequence. 

We can easily generalize from the bilayer case to mul- 






tilayer graphene with an arbitrary stacking order. Eigen¬ 
states at momentum orientation (f> k satisfy 

*fa) = (cue^.dBe^ 8 ^,-) 1 = Ufa)*(0), 

( 2 ) 

where 7*(0)tf(0) = e'F(O). Ufa) = Ufa)U{Q)U~ l fa) 
has matrix elements 

U ij fa)=Hi j ( 0)^“^, (3) 

and eigenvalues e that are independent of <p k . By com¬ 
paring the matrix elements in Eq. © with those in the 
original Hamiltonian, we can determine the phase fac¬ 
tor chirality parameters {Pi}- In general two sites con¬ 
nected by nearest-neighbor interlayer hopping have the 
same phase and within a layer Pb = Pa + 1- Using these 
two rules, {P;} is completely determined by the stacking 
sequence. Figure [T] illustrates their application in multi¬ 
layer structures with up to four layers. We explain below 
how these band structure properties influence electron- 
electron interaction physics in multilayer graphene sys- 


terns. 
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FIG. 1. Stacking diagrams and phase factor chirality pa¬ 
rameters {Pi} for (a) ABC and (b) ABA graphene, (c) 
Phase factors for all stacking arrangements from monolay¬ 
ers to tetralayers. We have chosen to set the phase factor of 
the sublattice 1A to zero. These results are for valley K. For 
valley K' the chirality parameters change sign. 

Exchange self-energy — Our goal in this paper is to ad¬ 
dress interaction effects in moderate carrier density mul¬ 
tilayer graphene systems, which are are weakly correlated 
two-dimensional Fermi liquids in which electron-electron 
interaction effects can be reliably addressed using pertur¬ 
bation theory. At leading order the electron self-energy 
is given by the unscreened exchange contribution: 

s ex (M) = -Y,J ( 4 ) 
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where f s , k is the Fermi function for the band s and 
wavevector fc, F kk , = |(s, k\s', k’)\ 2 is a wavefunction 
overlap factor, and V q = 2ne 2 /cq q is the two-dimensional 
Coulomb interaction. (We note that the Coulomb inter¬ 
action between layers with the layer separation d is given 
by V q {d) = V q e~ qd . Due to the small layer separation 
we approximate V q (d) ~ V q for the analytic calculations. 
The full numerical calculations using V q {d) do not change 
the results qualitatively.) It is conventional to absorb the 
self-energy at the Dirac point [k = 0) in the absence of 
carriers into the zero of energy. 

To understand the consequences for interaction physics 
of multilayer wavefunction chiral properties, it is instruc¬ 
tive to first consider the chiral two-dimensional elec¬ 
tron system (C2DES) Hamiltonians [S| that provide a 
low-energy effective theory of multilayer graphene. The 
Hamiltonian of a C2DES with the chirality index J is 



and yields eigenenergies e s , k = st± (hv\k\/t±) J , and 
eigenspinors |s, fe) = (s, e lJ< ^ fc ) t /\/2, where s = 
±1 for positive and negative energy states respec¬ 
tively. For a C2DES with the chirality J, F kk , = 
^ [1 + ss' cos J(4>k — = \ (1 + ss'n k ■ n k >), where 

n k = (cos J^fc,sin J<j> k ) is the pseudospin direction at 
k characterized by the chirality index J. Note that the 
overlap factor F kk , for a C2DES has the form of Heisen¬ 
berg interactions between pseudospins with orientation 
J(j) k and J(f> k '^§]- 

In Fig. H}a) intraband and interband contributions to 
the conduction band exchange self-energy of a C2DES are 
plotted. In Fig.[2}b) pseudospin chirality is illustrated by 
plotting the spin-1/2 pseudospin orientation appropriate 
for two-component spinors. As the chirality increases, 
the magnitude of each contribution is suppressed be¬ 
cause pseudospin orientation changes more rapidly with 
wavevector. Especially, the interband exchange is sup¬ 
pressed more strongly owing to the contribution from 
states occupying the infinite sea of negative energies. In 
Figs. m c) and (d) we compare C2DES exchange self¬ 
energies with those of ABC and ABA graphene multi¬ 
layers. At low carrier densities in ABC graphene, the 
relative chiral index of the dominant wavefunction com¬ 
ponents is 3 and the exchange self-energy resembles the 
weak form found in a C2DES with J = 3; as the density 
increases, interlayer hopping becomes less important, and 
the exchange self-energy eventually approaches that of a 
C2DES with J = 1. At low densities, ABA graphene 
is described Q by a direct product of chiral gases with 
J = 1 and J = 2. 

Density-density response functions and collective 
modes — Figure [3] plots loss functions Im[— e(q, w) -1 ] 
for several different multilayer graphene structures. Here 
e(g, to) is the dielectric function which we approximate 
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FIG. 2. (a) Exchange self-energies and (b) conduction band 

pseudospin direction (s = +1) for C2DESs with J = 1,2, 3, 4. 
Exchange self-energies of (c) the lowest conduction band 
(s = +1) of ABC trilayer graphene for n =j 10 11 ,10 13 cm -2 , 
and (d) the lowest (s = +1) and second lowest (s = +2) con¬ 
duction bands of ABA trilayer graphene for n = 10 12 cm -2 . 
(The bands of the ABA multilayer structures are shown in the 

inset.) Here £o = ——- and we use the effective fine structure 
2 

constant a = = 1 and momentum cutoff q c = 1/a. 


using the weak coupling random phase approximation 
(RPA) expression e(q,uj) = 1 — V q IIo(q,uj), and IIo(g, u>) 
is the non-interacting electron density-density response 
function: 


n 0 (q,w) = 9 sv ^2 

s,s' * 


^ k fs,k fs’.k+q ps,s' 

(2tt) 2 tiu + A*’ s fc+q + if] 


k,k-\-q ’ 


( 6 ) 

where g sv = g s g v = 4 is the spin-valley degeneracy, 
A fc fc +q = £s,k - £ s ',k+q, £ s ,k is the eigenenergy for the 
band index s and wavevector k , and i] is a positive in¬ 
finitesimal number. The black thick lines in Fig. [3] plot 
the boundaries of electron -hole continua within which 
ImIIo(q,w) is non-zero and electron hole excitations are 
allowed. When e(q, oj) = 0, the loss function has a 6- 
function peak corresponding to plasmon collective exci¬ 
tations. When the plasmon modes enter the electron 
hole continuum, they can decay into single electron hole 
pairs through the Landau damping process. In multi¬ 
layer graphene, plasmon modes decay through interband 
transitions. The shark-fin structures around ui = 0 re¬ 
flects from the chiral nature of the wavefunctions which 
lead to suppressed momentum-dependent scattering [7f. 

Ground State Energy — The ground-state energy is 
the sum of the non-interacting kinetic energy and inter¬ 
action (exchange-correlation) energies. The exchange- 
correlation energy can be expressed in terms of the 



FIG. 3. Loss function Im[— e(q,uj)~ 1 ] of (a) ABC, (b) ABA, 
(c) ABCA, and (d) ABAB stacked multilayer graphene for 
n = 10 12 cm -2 and a = 1 with y = 5 x 10 -5 £f. The thick 
black lines indicate boundaries of the electron-hole continua 
and the insets in each panel show the energy band struc¬ 
ture. In the ABA structure and in other multilayer struc¬ 
tures with mirror symmetry, some interband transitions do 
not contribute to plasmon Landau damping, as indicated by 
dotted lines in panel (b). 


density-density response function^ by applying the 
integration-over-coupling-constant method and appeal¬ 
ing to the fluctuation-dissipation theorem. The RPA ap¬ 
proximation to the exchange-correlation energy is justi¬ 
fied in part by the relatively large spin-valley flavor de¬ 
generacy g sv = 4 which makes the RPA bubble-diagram 
contributions to the energy more dominant Q. For tech¬ 
nical reasons it is convenient to separate the first-order 
exchange-correction to the interaction energy and higher 
order corrections commonly referred to as the corre¬ 
lation energy. The dependence of the exchange and 
RPA correlation energies on carrier density can then 
be expressed^), |Io| as integrals along the imaginary fre¬ 
quency axis: 


£ex — 


£corr — 


R f £q r°°dv , . . 

i f d 2 q 


( 7 ) 


+ In 


0 dv 

2 n J (27t) 2 J 0 n 
1 - V q n 0 (q,iu) 


V q 5U.o{q,iv) 


1 - V q n 0 (g,ii/)| n=0 


where SHo(q,w) = IIo (q,w) — IIo(q, iv)\ n — 0 . We use 
the momentum cutoff q c = 1 /a to remove the ultraviolet 
divergences in the momentum integrals. 

Using these expressions, we find that in terms of the di¬ 
mensionless coupling constant ap = e 2 /eoh v F = (v/vp)a 
where Vf is the Fermi velocity, the exchange energy is 
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given by 

e 2 

e ex = hvFkpCicXF = — kpC±, (8) 

£o 

and the correlation energy for small «p has the form of 
Ecorr = kvpkF (C^QJp + ••■), (9) 

whereas in the strong coupling limit (ap 1), 

Ecorr = Hv F k F (HiCKf + H 0 H-) . (10) 

The coefficients {Ci} and {H,} in these expressions have 
weak density dependence through the response function 
and the momentum cutoff. 



FIG. 4. Comparison between the ground-state exchange 
energies (left panel) and correlation energies (right panel) of 
several different multilayer graphene (thick lines) structures 
and C2DES systems (dotted lines) as a function of carrier 
density for a = 0.05. 

To address multilayers, we first discuss the case of 
C2DES models whose exchange and correlation energies 
exhibit a systematic dependence on their chirality index 
J. As we see from Eq. m that the exchange energy is 
approximately proportional to kp irrespective of the in¬ 
teraction strength. The positive value of C\ for J = 1 
reflects the dominance in this case of interband exchange. 
For J > 1, the interband contribution is suppressed be¬ 
cause of the larger chirality indices and the exchange en¬ 
ergy turns negative (C\ < 0). For the correlation en¬ 
ergy in the strong coupling limit, we find from Eq. USD 

that £ CO rr = (e 2 /eo)fcF-Di + JepD 0 -|-, where £f is the 

Fermi energy. Note that C\ = —Hi; thus for J > 1 £ corr 
is positive (Hi > 0), whereas for J = 1 £ corr is negative 
(Hi < 0). In the weak coupling limit, we see from Eq. © 
that £ cor r oc t^kpap k^' 1 - 

Figured] illustrates these properties of C2DES models 
and compares these exchange-correlation energy proper¬ 
ties with those of multilayer graphene systems. For multi¬ 
layer graphene, at low densities, the exchange and corre¬ 
lation energies follow those of the largest J C2DES model 


contained within its low-energy bands because they are 
responsible for the largest density of states (DOS). For 
example, for ABA stacking the exchange and correlation 
energies at low densities follow those of a J = 2 C2DES 
because a J = 2 C2DES has a larger DOS than a J = 1 
C2DES. As the carrier density increases, interlayer hop¬ 
ping becomes less important and the exchange and cor¬ 
relation energies begin to approach those of monolayer 
graphene. 

Summary and Discussion — In this paper we have ex¬ 
ploited the simple dependence of band wavefunctions on 
momentum orientation to simplify many-electron pertur¬ 
bation theory calculations for multilayer graphene, and to 
bring out the relationship between stacking sequence and 
carrier-carrier interaction phenomena in this interesting 
class of materials. By explicit calculations for a variety of 
different structures we have shown that the exchange self¬ 
energies and related exchange-correlation energy features 
in multilayer graphene systems follow those of C2DES 
models at low carrier densities, but cross over to be more 
similar to those of monolayer graphene as carrier densi¬ 
ties increase. The rotational transformation of the chiral 
wavefunction is very general and can be applied even in 
the presence of site energy variations or remote hopping 
terms, unless momentum-dependent hopping terms do 
not appear in the Hamiltonian. For example, if the re¬ 
mote interlayer hopping term 72 is included in the Hamil¬ 
tonian, it modifies {c^}, but not the angular part of the 
wavefunction in Eq. ©■ 

The model we employ, however, does not include 
momentum-dependent remote interlayer hopping terms. 
We also use the weak-coupling RPA in the calculation. 
Both limit the applicability of our calculations to mod¬ 
erate to high carrier densities, and at very low carrier 
densities correlations frequently become strong and lead 
to broken symmetry ground states not captured by the 
RPA [13 . 

Our theory, however, captures important observable 
effects produced by interactions which strongly depend 
on the stacking sequences such as plasmon collective ex¬ 
citations and self-energies. The dependence of ground 
state energies on carrier and spin densities are responsi¬ 
ble for renormalized electronic compressibility and spin 
susceptibility, respectively. For example, the electronic 
compressibility n is given by k -1 = n 2 dfi/dn , where 
/r = d{netot)/dn is the chemical potential of the interact¬ 
ing system, and £t 0 t is the total ground-state energy per 
particle. In an ordinary parabolic band two-dimensional 
electron system the compressibility famously becomes 
negative at low carrier densities [ 12y . In contrast, for a 
C2DES we find that in the low-density strong-coupling 
limit, kq/k 1 + J(J + 2)Do/2 with Ho > 0 for J = 1 
and Ho < 0 for J > 1. (Here Kg is the non-interacting 
compressibility.) It follows that for a J = 1 C2DES, the 
electronic compressibility is strongly suppressed by inter¬ 
actions due to the interband exchange contribution and 
remains positivepjl. For a J > 1 C2DES, the interband 
exchange contribution is suppressed with the chirality; 
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thus, the electronic compressibility is enhanced by the 
interactions. Interestingly, we find that for J > 5, the 
compressibility can be negative in the low-density limit, 
suggesting an instability toward other ground states, 
whereas in the absence of correlation, this occurs for 
J > 2. For multilayer graphene, at low carrier densities, 
the compressibility follows the trend of the correspond¬ 
ing C2DES, but as the carrier density increases, the com¬ 
pressibility follows that of monolayer graphene with sup¬ 
pressed compressibility [H, H3 |. showing non-monotonic 
behavior arising from competition between the intraband 
exchange, interband exchange, and correlation. 

In conclusion, our new approach allows us to effectively 
calculate the quasiparticle and thermodynamic proper¬ 
ties of interacting many-body chiral systems. We show 
that as the chirality increases, the exchange contribu¬ 
tion to the single particle energy is suppressed and the 
correlation contribution increases, indicating that the 
exchange-correlation is controlled by the stacking ar¬ 


rangement. Our results suggest that correlation effects 
play a more important role in a system with a large chi¬ 
rality; thus, we expect that rhombohedral graphene with 
periodic ABC stacking could show exotic interaction- 
induced phenomena such as ordered states and non-Fermi 
liquid behavior (Tlj| . 

ACKNOWLEDGMENTS 

This research was supported by the Basic Science Re¬ 
search Program through the National Research Foun¬ 
dation of Korea (NRF) funded by the Ministry of 
Science, ICT and Future Planning under Grant No. 
2012R1A1A1013963 (YJ and HM), Basic Science Re¬ 
search Program 2009-0083540 (EHH) and DOE Materi¬ 
als Science and Engineering grant DE-FG03-02ER45958 
and by Welch Foundation grant TBF1473 (AHM). HM 
thanks Eun-Gook Moon for helpful discussions. 


[1] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. 
Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 
(2009). 

[2] S. Das Sarma, S. Adam, E. H. Hwang, and E. Rossi, Rev. 
Mod. Phys. 83, 407 (2011). 

[3] D.N. Basov, M.M. Fogler, A. Lanzara, F. Wang, and Y. 
Zhang, Rev. Mod. Phys. 86, 959 (2014). 

[4] H. Raza (Ed.), Graphene Nanoelectronics, Springer 

( 2012 ). 

[5] H. Min and A. H. MacDonald, Phys. Rev. B 77, 155416 
(2008); H. Min and A. H. MacDonald, Prog. Tlieor. Phys. 
Suppl. 176, 227 (2008). 

[6] H. Min, G. Borghi, M. Polini, and A. H. MacDonald, 
Phys. Rev. B 77, 041407(R) (2008). 

[7] G. Borghi, M. Polini, R. Asgari, and A. H. MacDonald, 
Phys. Rev. B 80, 241402(R) (2009). 


[8] G. F. Giuliani and G. Vignale, Quantum theory of the 
electron liquid, Cambridge University Press (2005). 

[9] M. Polini, R. Asgari, Y. Barlas, T. Pereg-Barnea, and A. 
H. MacDonald, Solid State Commun. 143, 58 (2007). 

[10] Y. Barlas, T. Pereg-Barnea, M. Polini, R. Asgari, and A. 
H. MacDonald, Phys. Rev. Lett. 98, 236601 (2007). 

[11] V. N. Kotov, B. Uchoa, V. M. Pereira, F. Guinea, and 
A. H. Castro Neto, Rev. Mod. Phys. 84, 1067 (2012). 

[12] J. P. Eisenstein, L. N. Pfeiffer, and K. W. West, Phys. 
Rev. Lett. 68, 674 (1992). 

[13] S. Viola Kusminskiy, J. Nilsson, D. K. Campbell, and A. 
H. Castro Neto, Phys. Rev. Lett. 100, 106805 (2008). 

[14] G. Borghi, M. Polini, R. Asgari, and A. H. MacDonald, 
Phys. Rev. B 82, 155403 (2010). 



